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ABSTRACT 

Galaxy centers are residing places for Super Massive Black Holes (SMBHs). Galaxy mergers bring 
SMBHs close together to form gravitationally bound binary systems which, if able to coalesce in 
less than a Hubble time, would be one of the most promising sources of gravitational waves for the 
Laser Interferometer Space Antenna (LISA). In spherical galaxy models, SMBH binaries stall at a 
separation of approximately one parsec, leading to the "final parsec problem" (FPP). On the other 
hand, it has been shown that merger-induced triaxiality of the remnant in equal-mass mergers is 
capable of supporting a constant supply of stars on so-called centrophilic orbits that interact with the 
binary and thus avoid the FPP. In this paper, using a set of direct A^-body simulations of mergers of 
initially spherically symmetric galaxies with different mass ratios, we show that the merger-induced 
triaxiality is also able to drive unequal-mass SMBH binaries to coalescence. The binary hardening 
rates are high and depend only weakly on the mass ratios of SMBHs for a wide range of mass ratios q. 
There is, however, an abrupt transition in the hardening rates for mergers with mass ratios somewhere 
between q ~ 0.05 and 0.1, resulting from the monotonic decrease of merger-induced triaxiality with 
mass ratio q, as the secondary galaxy becomes too small and light to significantly perturb the primary, 
i.e., the more massive one. The hardening rates are significantly higher for galaxies having steep cusps 
in comparison with those having shallow cups at centers. The evolution of the binary SMBH leads to 
relatively shallower inner slopes at the centers of the merger remnants. The stellar mass displaced by 
the SMBH binary on its way to coalescence is ~ 1 — 5 times the combined mass of binary SMBHs. The 
coalescence time scales for SMBH binary with mass ~ lO^M© are less than 1 Gyr and for those at the 
upper end of SMBH masses 10^ Mq are 1-2 Gyr for less eccentric binaries whereas less than 1 Gyr for 
highly eccentric binaries. SMBH binaries are thus expected to be promising sources of gravitational 
waves at low and high redshifts. 

Subject headings: Stellar dynamics - black hole physics - Galaxies: kinematics and dynamics - Galaxy: 
center. 



1. INTRODUCTION 

Super Massive Black Holes (SMBHs) are now a well 
established component o f galaxies having a sizable bulge 
(jFerrarese fc Fo rd"2005l. In the paradigm of the A Cold 
Dark Matter cosmology, galaxies are formed via hierar- 
chical merging. If more than one of these merging galax- 
ies contain a SMBH, the formation of a b inary SMBH 
is almost inevitable (jBegelman eFall [19801) . There are 
cases in which there is a clear observational evidence for 
two widely separated SMBHs as well as s ome circumstan - 
tial evidence for a true SMBH binary (jKomoss al [2006h . 
Upon coalescence, SMBH binaries would be the highest 
signal-to-noise ratio sources of gravitational waves, de- 
tectable from essentially any cosmological redshift, for 
future space-borne gravitational wave detectors such as 
the Laser Interferometer Space Antenna (L ISA) ( Hughesl 
l2003HBarack fc Cutled[2004l: iDanzmann et al. 201% 

The paradigm for SMBH binary evolution, after a 
merger of gas-poor galaxies, consists of three distinct 
phases (jBegelman et all I198O0 . First, the two SMBHs 



sink towards the center due to the dynamical friction ex- 
erted by the stars and the dark matter. Dynamical fric- 
tion acts together with "gravitational slingshot effect" as 
SMBHs form a bound pair with semi- major axis a ^ r^, 
where rh is the binary's infiuence radius, which typically 
is defined to be the radius which encloses twice the mass 
of the binary in stars M(< rh) = 2M,, where M, is the 
mass of the binary. Dynamical friction stops being an 
effective driver of inspiral whe n the binary reaches th e 
hard binary separation a Oh (|Quinlanlll996l: 1^312001 



where fi^ is the binary's reduced mass and a is the ID 
velocity dispersion. In the second phase, slingshot ejec- 
tion of stars is the dominant mechanism for taking en- 
ergy and angular momentum away from SMBH binary. 
Thirdly, the binary eventually reaches a separation acw 
at which the loss of orbital energy due to Gravitational 
Wave (GW) emission drives the final coalescence. The 
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transition from the first to the second phase is prompt 
provided that t he mass ratio of the r emnants is not too 
smah q ^ 0.1 (iCahegari et aLir2011[ ). In contrast, the 
subsequent transition from the second to the third phase 
could constitute a bottleneck for the binary evolution 
towards final coalescence caused by lack of stars which 
can interact with the binary. This is the so-called Final 
P arsec Pro blem. 

lYul (|2002[ ) early pointed out, based on an analysis based 
on the HST sample of ne arby elliptical galaxies and spiral 
bulges from lFaber et alJ ( jT997l ), that fiattening and non- 
axisymmetries would help stella r dynamics in b r inging 
SMBH binaries to coalescence. iMerritt fc PoonI (|20o4 ) 
built self-consistent cuspy, triaxial models with a single 
SMBH at the center and showed that such models could 
be concocted with a significant fraction of centrophilic or- 
bits that would efficiently drive the hardening rate of a bi- 
nary if it were present at the center. Berczik et al.. (.2006,) 
showed, for the first time with iV-body simulations of ro- 
tating King models, that the FPP could potentially be 
solved by purely stellar dynamical models that develop a 
significant amount of triaxiality. By including the post- 
Newtonian CPAQ_termsJntotlK of motion of 
the binary, iBerentzen et al.l (|2009( ) explicitly showed that 
the coalescence times, if the latter models are scaled to 
real galaxie s , are clearly shorter than a Hubble time. 
iPreto et al.l ()2009[ ) presented preliminary results sugges- 
tive that merger-induced triaxiality resulting from the 
merger of spherically symmetric models would lead to 
qualitatively similar results - namel y the prompt coa - 
lescence of resulting S MBH binarie s . |Khan et al.l ()2011[ ) 
(hereafter paper I) and IPreto et al.l ()2011[ ) (hereafter pa- 
per H) confirmed and generalized the latter results by 
showing that galaxies which form via equal mass merg- 
ers are mildly triaxial. They have shown that hardening 
rates of SMBH binaries resulting from equal-mass merg- 
ers are substantially higher than those found in spher- 
ical nuclei, and are independent of the total number 
N of stars, thus allowing them to extrapolate results 
to real galaxies. Presumably, the nonspherical shapes 
of merger remnants reported in these studies result in 
a large population of stars on centrophilic orbits. In- 
side the influence radius of SMBH, the centrophilic orbit 
family incl udes saucer or cone orbit s in the axisymmetric 
potentials (jSridhar fc Toumalll999[ ) and pyramids orbits 
(jugular momentum close to zero) in triaxial potentials 
([Mcrritt & Vasiliev 2011). Outside the infiuence radius, 
most of t he centrophilic orbits a re choatic in a triaxial 
potential (| Poon fc Merrittil2001[ ). 

The inspiral of a binary SMBH is expected to leave a 
characteristic imprint in the morphological and dynam- 
ical properties of the newly formed galactic nucleus fol- 
lowing t he merger: e.g., th e bending and precession of ra- 
dio jets (jRoos et al.lll993D, the variabil ity in optical light 
curve of quasar (iValtonen et al.ll2008[ ). X-shaped radio 
lobes ()Merritt fc Eker£ll2002D and, with some likelihood, 
a partially deplete d core instead of a steep cusp in Ijright 
elliptical galaxies (|Grahamll20"ol iKormendv et al.ll2009| ) 
could be explained by different mod els of SMBH binary 
evolution prior or after coa lescence (jMerritt et al.l 120071 : 
iGualandris fc Merrittil2008l ). 

Detailed surface photometry of ellipticals has revealed 
that their surface brightness is well described by a Sersic 



profile, log/(r) oc r^/", over most of the body of a 
bulge or early-type galaxy. However, systematic devia- 
tions from this profile are found in the innermost regions 
close to the central SMBH - either in the form of a light 
(mass) excess or deficit with respect to that would result 
from t he extrapolation of the Sersic profile towards the 
center (|Kormendv et al.ll2009[ ). In other words, observa- 
tions seem to reveal a dichotomy in the central density 
profiles of bulges and early-type galaxies: giant, high- 
luminosity objects have relatively shallow central density 
profiles while normal, low-luminosity ones show steeper 
density profiles in their center. The former are often 
called core galaxies and show density profiles with central 
logarithmic slopes 7 < 1, and typically harbor SMBH 
of mass > 10^ Mq; while the latter are normally called 
power law galaxies, have central 7 > 1 and lighter SMBH 
< few X 1O^M0. The explanation of cores is non-trivial 
since dry mergers should preserve the highest density in 
phase space given that, by virtue of being devoid of ap- 
preciable a mounts of gas, their dynam ics is essentially 
collisionless (jBinnev fc Tremainell2008[ ). 

Galaxies having a single mass stellar population are 
expected to have Bahcall fc Wolf (1976 ) cusp around 
the SMBH with p{r) oc r-",a ~ 7/4 for r < rh af- 
ter a relaxtion time (|Preto et al.l I2004D . In the more 
realistic case where there is a range of stellar masses, 
less massive objects follow a shallow profile, a ~ 3/2, 
while heavier objects (e.g. stellar blac k holes) develop 
steeper cusps a ^ 2 around the S MBH (Bahcall fc Wolj 
Il977t I Alexander fc Hopman I [20091 Whether a given, 
or the typical, nucleus has reached this relaxed state 
or not depends on its "initial" profile and on how long 
ago it was formed. Different studies have reached differ- 
ent conclusions regarding the most probable state of the 
stellar distribution around the S MBH after a relaxation 
time. IPreto fc Amaro-Seoand ()2010f ). who performed 
Fokker-Planck and A''— body simulations of a two com- 
ponent galaxy nucleus model made of solar mass stars 
and stellar mass black holes, obtained mass-segregated 
stellar cusps in ^ 1 / 4 of a relaxation time. However, 
IGualandris fc MerrittI (|2012| ) state roughly the opposite 
conclusion: that even after a few central relaxation times, 
the distribution of stars can be very different than in the 
steady-state, mass-segregated solutions. 

The analysis of the number counts of spectroscopically 
identified, old sta rs in the s ub-parscc region of o ur own 
Milky Way (Buch holz etrall !2009: D o et al.ll200"l seems 
to favor the presence of a very shallow profile for the vis- 
ible fraction of the old stellar population. This is consis- 
tent with the finding that the two-body relaxation time 
at the influence radius of the SMBH in the center of the 
Milky Way is longer than a H ubble time (Mcrritt 201(| 
IPreto fc Amaro-Seoand[2010f l. On the other hand, low- 
luminosity spheroids often contain nuclear star clusters 
which ha ve central rela xation times shorter than a Hub- 
ble time (|Merrittll2009t) . 

In dry mergers, i.e. in the absence of a dynamically 
significant amount of gas, the slingshot ejection of stars 
by the central SMBH binary drives the binary inspiral. 
It has early been proposed that inspiraling binaries could 
carve a core inside ~ rh, thus providing an explanation 
for the mass d eficits present in the center of gas-poor gi- 
ant ellipticals (|Makino fc Ebisuzakilll996l) . On energetic 
grounds it can be shown that in order for the binary to 
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reach coalescence, it needs to eject an amount of stel- 
lar mass of order of its own mas s, Mei ^ a x M«, where 
a = 0(l) and M, = M,^ +M.2 (iMerritt fc Milosa vlievid 
[2005t iPerets k Alexanded l2008nnMCTrirt (2006a) stud- 
ied the mass deficits created by SMBH inspirals evolv- 
ing in spherical symmetric nuclei and found an average 
a ^ 0.5. However, unless their mass M, < few x 10^ Mq, 
the evolution of SMBH binaries in gas-poor spherical 
galaxies tends to stall due to depletion, on the (local) 
dynamical timcscale, of the pool of stars whose orbits 



intersect the binary (e. g. [Milosavlievi c & Merritt 2003; 
Makino fc Fun ato 2004': lBerezik et all 120051: IKhan et al] 
201lHPreto'eral...2011.) . The binaries studied bv lMerrittj 



( 2006afl ' did not reach coalescence due to the FPP. There- 
fore they estimated mass deficits at the time when a hard 
binary forms in their simulations. iMerritt et al.l ()2007l ) 
extended the latter study by following the inspiral of bi- 
naries up to coalescence in the lower mass range using an 
approximate Fokker-Planck description for stellar scat- 
tering. They find substantially higher mass deficits than 
before, albeit restricting their calculations to spherical 
models and very shallow cusps (with initial 7 ~ 1/2) in 
strong contrast with the larger 7 typical of compact low- 
mass nuclei. More precise estimates of a from more real- 
istic galaxy merger studies are therefore certainly needed; 
in particular, ones that do not rely on any of the follow- 
ing approximations: spherical symmetry, low 7, crude 
model for inferring the mass ejected from the increase in 
the binary's orbital energy, or those generally inherent 
to the Fokker-Planck formalism. In this respect, prelim- 
inary studies of mass deficits created by the merger of 
equal mass galaxies, with moderate SMBH mass ratio, 
moderate central slope 7 ~ 1, and that generate mildly 
triaxial remnants, suggest that the resulting mass deficits 
could be substant ially higher Mcj ^ 2M, within < 3rh 
(jPreto et al.ll2009l ) - a result which we confirm and gen- 
eralize in this paper. Note, however, that alternative 
scenarios for explaining the formation of cores have been 
proposed (|Gualandris fc Merrittil2008[ ). 

In this paper, we revisit these problems by detailed 
iV-body simulations using a range of more realistic mod- 
els of merging galactic nuclei with varying mass ratio q 
and central logarithmic slopes 7. We restrict the simu- 
lations to g = Mgai,s/-/Wgai,p — A/,2/Af,i. We study the 
shape of the merger remnant of galaxies having different 
density profiles and different mass ratios. For each of 
our models, the hardening rate of the SMBH binary, the 
axis ratios of the merger remnant and the resulting mass 
deficit are calculated. In section 2, we describe our initial 
models and numerical methods used in our study. The 
evolution of SMBH binary in galaxy mergers is described 
in section 3, as well as the time evolving densities and 
shapes of newly formed galaxies resulting from the galaxy 
merger. Section 4 describes the analysis and results of 
"mass-deficits" resulting from the galaxy mergers. The 
coalescence timescales for binary SMBHs from the time 
of formation of binary till the full coalescence of SMBHs 
due to emission of gravitational waves are presented in 
section 5. Finally, section 6 summarizes and discusses 
the results of the paper. 

2. INITIAL CONDITIONS AND NUMERICAL METHODS 
2.0.1. The host galaxies and their SMBHs 



Our aim is to study the dependence of the harden- 
ing rates with the mass ratios and central concentra- 
tion of the nucleus - as well as the resulting imprint of 
the SMBH inspiral on the stellar distributions of post- 
coalescence nucleus. Following closely the initial set-up 
of our previous numerical experiments (see papers I and 
II), we represent the individual galaxies or galactic nuclei 
by spherically symmetric TV-body models. The density 
profile of the models follows the Dehnen (1993) profile of 
the form 



(3 - 7)Mgal 



4tt 



(2) 



where Mgai is the total mass of the galaxy, ro is the 
scale radius, and 7 is the inner logarithmic slope. For 
r rp, the density declines with a logarithmic slope of 
—4, which, although not close to the exterior density of 
real galaxies, provides a convenient cut-off that avoids 
the waste of CPU time in following stars whose orbits 
are too far from the SMBH binary to have any signif- 
icant impact on its long term evolution. The resulting 
cumulative mass profile is given by 



MD(r) = Mg 



3-7 



(3) 



In addition, we represent the SMBHs by massive parti- 
cles with zero velocity placed at the center of both the 
primary and the secondary galaxies. The masses of the 
SMBHs are set to be 0.5 percent of their host galaxy. 
This value is few times higher th an the currently ac- 
cepted value o f M^^i/M, ^ 0.001 ([F errarese fc Merritd 
2000; Gcbhar dt et al.l 120001: IMerritt fc Ferrarese 200ll 
HiirinjjLRa 12004'). However, some galaxies may have 
larger value of this ratio than others. 

Thus the two SMBHs in the simulations have the same 
mass ratio q as their host galaxies. The Dehnen models 
have an isotropic distribution function of velocities (no 
net rotation) , are spherically symmetric and are initially 
set-up to be in dynamical equilibrium with a gravita- 
tional potential given by <I>(r) = — GM,/r-|-$*(r), where 
<i>*(r) is the gravitational potential due to the stars alone. 

The galaxies size ratio scales with the corresponding 
mass ratio as R2/R1 oc y/ M2/M1 . By choosing four 
different values of 7 = 0.5, 1.0, 1.5, 1.75, the central part 
of our galaxy models represent the variety of observed 
density profiles in both early and late type galaxies. 

For Dehnen models the ratio of the effective radius 
(projected half-mass radius) Rc to the half-mass radius 
fi/2 is given approximately as 



Re 

ri/2 



0.75, 



(4) 



and shows only a weak dependence on 7. In the follow- 
ing we refer to the more massive galaxy as the primary 
galaxy and the lighter one as the secondary. In our model 
units we use for our primary galaxy a total mass and scale 
radius of Mgai.p = ro,p = 1. We also set the gravitational 
constant to G = 1. 

In the absence of relativistic effects on the binary's 
orbit, our models are scale-free and can be a posteriori 
scaled to galaxies of different total mass and size. How- 
ever, the inclusion of radiation reaction effects on the 
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TABLE 1 

Model parameter and properties of the primary galaxy. 



Model 


N 


7 


ri/2 


Reft 


A 


128k 


0.5 


3.13 


2.35 


B 


128k 


1.0 


2.41 


1.81 


C 


128k 


1.5 


1.69 


1.27 


D 


64k 


1.75 


1.35 


1.01 



Note. — Columns from left to right show the model name, 
the particle number A'^, the half-mass radius ri^2 s-iid the effective 
radius ReS, respectively, for our primary galaxy models. 

binary due to GW emission introduces an absolute scale 
associated to the universal value of the speed of light. 
Such scalings are described in detail in section 5 when 
coalescence times are explicitly computed. 

The effective radii of our models can be compared to 
the effective radii in observed galaxies to calculate tq. 
By fixing galaxy mass Mgai our models can be scaled to 
physical units using following relations: 



/ ]\ T \ 1/2 / \ -1/2 

We have performed a series of A^-body experiments 
where we vary the mass ratio of the galaxies and SMBHs 
q and the inner density slope 7. In Table [T] we sum- 
marize the model parameters and properties for the pri- 
mary galaxies used in our simulations. The computa- 
tional wallclock time increases with increasing cuspiness 
(higher 7) of the density profile since particles in the 
center need very small numerical timesteps in order to 
resolve accurately the (locally) strongly varying gradi- 
ents of the gravitational potential. In order to keep the 
computational time within reasonable limits we built our 
model D (highest 7 = 7/4) with only half as many par- 
ticles as the other models. 

The secondary galaxies used in our merger simulations 
have the same density profile, i.e., the same 7, as the pri- 
mary galaxies, but have different masses Mgai,s and scale 
radii tq^s- Both primary and secondary galaxy models 
have equal number of particles N. 

The initial center of mass positions and velocities for 
the two galaxies are calculated from the Keplerian orbit 
of the equivalent two-body problem, with given apo- and 
peri-centers, Ta and rp, respectively. The two galaxies 
start at the apo-center of their relative orbit. The initial 
separation between the center of mass of the two galaxy 
is Ar = 15ro,p; and since the half-mass radius of each 
nucleus is i?i/2 < 2.5, this choice ensures that the galax- 
ies are initially well separated while, at the same time, 
we minimize computing time. The initial relative veloc- 
ity of the two galaxies is chosen such that the SMBH 
separation at first peri-center passage is ~ 2ro,p; in other 



TABLE 2 

Parameters of the galaxy mergers study 



Run 


Aftot 




n 
H 


' U,s 


Al 


256k 


0.5 


0.1 


0.316 


A2 


256k 


0.5 


0.25 


0.5 


A3 


256k 


0.5 


0.5 


0.707 


A4 


256k 


0.5 


1.0 


1.0 


Bl 


256k 


1.0 


0.1 


0.316 


B2 


256k 


1.0 


0.25 


0.5 


B3 


256k 


1.0 


0.5 


0.707 


B4 


256k 


1 n 


1 n 


1 n 


B5 


256k 


1.0 


0.05 


0.22 


CI 


256k 


1.5 


0.1 


0.316 


C2 


256k 


1.5 


0.25 


0.5 


C3 


256k 


1.5 


0.5 


0.707 


C4 


256k 


1.5 


1.0 


1.0 


Dl 


128k 


1.75 


0.1 


0.316 


D2 


128k 


1.75 


0.25 


0.5 


D3 


128k 


1.75 


0.5 


0.707 


D4 


128k 


1.75 


1.0 


1.0 



Note. — Col.(l) Galaxy merger model. Col(2) Number of par- 
ticles. Col(3) Central density profile index 7. Col(4) Mass ratio 
of merging galaxies and their black holes. Col(5) Scale radius of 
secondary galaxy. 

words, the initial orbit of the binary galaxy has eccentric- 
ity 0.75, which corresponds to a circularity parameter 
value of e = L/Lc 0.66 in rough agreement with the 
values measured from f ull galaxy scale merger simula- 
tions (|Kazantzidisl[20l0l ). 

2.1. Numerical Methods and Hardware 

The iV-body integrations are carried out using 
the 0GRAPEQ (jHarfst et all I2007D . a parallel, direct- 
summation iV-body code that uses general purpose hard- 
ware, namely Graphic Processing Units (GPU) cards 
to accelerate the computation of pairwise gravitational 
forces between all particles. Our updated version of 
the (/)GRAPE (— Parallel Hermite Integration with 
GRAPE) code uses the 4-th order Hermite integration 
scheme to advance all particles with hierarchical individ- 
ual block time-steps, together with the parallel usage of 
GPU cards to calculate the acceleration a and its first 
time derivative a between all p articles. For the fo rce 
calculation we use the SAPPORO (|Gaburov et al.ll2009D li- 
brary which emulates on the GPU the standard GRAPE- 
6 library calls. 

As (/)GRAPE does n ot include the regularization 
(|Mikkola fc AarsethI I1998D of close encounters or bina- 
ries, we have to use a gravitational (Plummer-)softening 
between all components. The softening length is cho- 
sen to be sufficiently small (i.e. 5 x 10~^ in our model 
units) to keep the (dense) stellar system as collisional as 
possible/necessary. In the few runs in which we include 
the VAf terms to the equations of motion of the SMBH 
binary, the softening between the SMBH particles is set 
equal to zero. Our impl ementation of th e VAf terms is 
identical to that used in IBerentzen et al.l (j2009D , except 
that we cut-off the VN" expansion at the next order. We 
stress that the binary is, in all cases, integrated taking 
into account the gravitational field from the stellar clus- 
ter. 

The iV-body integrations were carried out on three 

^ |f tp : //ftp ■ arl ■ unl-heldelberg ■ de/staf f /berczlk/phl-GRAPE/ 1 
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computer clusters. "titan", at the Astronomisches 
Rechen-Institut in Heidelberg, the GPU-enhanced clus- 
ter "kolob" at the University of Heidelberg and "laohu" 
in Beijing. 

For a detailed discussion on the iV-body methods and 
choice of parameters, we refer the reader to paper I. 

3. SMBH BINARY EVOLUTION IN GALAXY MERGERS 

Figure [1] shows the relative separation R between the 
two black holes during a galaxy merger. According to the 
equivalent two-body trajectory used to set up the initial 
galaxy orbital parameters, the galaxy centers, and hence 
the SMBHs, should reach a separation of 2ro,p during the 
first pericenter passage. However, as figure [T] shows, in 
case q > 1/2, the SMBH separation shrinks well below 
2ro,p during the galaxy merger phase. As q decreases, 
and hence the mass of the secondary galaxy decreases as 
well, the galaxies reach the first pericenter passage after 
following more closely the corresponding two-body orbit. 
In the same figure the color of the arrows signals the time 
T at which the two SMBHs in a particular model form 
a binary system. As q decreases, the arrows move from 
left to right signaling the longer interval of time between 
the start of merger and formation of the binary. Neither 
the time T nor i?pori are as sensitive to 7 as they are to 

q- 

When the SMBHs become bound their separation is 
i? ~ Th. Once a SMBH binary is formed, dynamical fric- 
tion and the ejection of stars by the gravitational sling- 
shot, act together to efficiently extract angular momen- 
tum from the massive binary and the binary separation 
shrinks rapidly. When the semi-major axis a « ah, this 
initial rapid phase of binary hardening comes to an end 
as the pool of stars whose orbits cross the binary orbit 
gets depleted on the (local) dynamical timescaleo 

The subsequent hardening of the binary occurs only 
if the loss cone orbits are replenished^ In a spherical 
galaxy with no gas, the only dynamical mechanism for 
replenishing loss cone orbits is two-body relaxation, with 
the corresponding timescale for re populating the loss 
cone orbits scaling as ~ N/ In N (jBinnev fc Tremaind 
[2008I) . and the inspiral rate thus becomes strongly N - 
dependent (|Makino fc Funatoll2004l : IBerczik et al.ll2005D . 
In papers I and II, it was shown that the hardening phase, 
in the case of a mildly triaxial remnant nuclei result- 
ing from a galaxy merger, is essentially iV-independent 
- which allows our results to be scaled to real galax- 
ies whose typical number of stars, iVstar, is always sev- 
eral orders of magnitude larger than the maximum num- 
ber reachable with state-of-the-art direct TV-body simu- 
lations. Furthermore, it was also shown in Papers I & 
II that the flux of stars into the loss cone in a triaxial 
remnant is also consistently higher than in the spherical 
case, meaning the binary evolution and coalescence can 
occur much faster than in spherical models. 

Figure [2] shows the time evolution of the binary's in- 
verse semi- major axis 1/a from the time T when the two 

^ The semi-major axis a and eccentricity e of the binary are 
defined here via the standard Keplerian relations, i.e., neglecting 
effects of the field stars. 

The loss cone is the region of phase space corresponding, 
roughly speaking, to orbits that cross the b inary, i.e. with an- 
gular m omentu m J < Jj^ = ^/(JW^Ja^~, where / = 0(1) 
IlLightman fc ShapirollmW) . 



black holes become bound. An initial phase with faster 
binary hardening becomes more evident for steep inner 
density profiles (7 = 1.5, 1.75) - presumably correspond- 
ing to the clearing of the original loss cone (Yu 2002). 
The binary hardening rates s are calculated by fitting 
straight lines to a^^{t) in the late phase of binary evolu- 
tion from T = 60 — 100. The results are shown in Figure 
[3l The hardening rates depend very weakly on the mass 
ratio of binary SMBHs, in contrast with the ~ M, de- 
pendence found in equal mass merger study (papers I & 
II) and, for fixed Af,, on the mass ratio q. On the other 
hand, the value of s{t) increases significantly for higher 
values of 7. Note that for the single run with the smallest 
mass ratio in our sample, q = 0.05 and 7 = 1, the hard- 
ening rates are significantly weaker than the remaining 
7=1 models. 

Following Paper II, the SMBH binary hardening rate, 
which is a measure for the energy loss of the binary, is 
given by 

where J-\c{E,t) is stellar flux into the loss cone, (C) = 
0(1) is a dimensionless quan tity obtained f rom three- 
body scattering experiments (jOuinlanl 119961 ) and E = 
GM,/r -I- $*(r) — 1/2 is the (specific) energy of each 
star {E > for bound stars) . The behavior of the hard- 
ening rates can thus be qualitatively understood as fol- 
lows. The flux Tic{E,t), and its time evolution, depends 
on the degree of symmetry in the underlying gravita- 
tional potential (or, equivalently, the orbit families sup- 
ported by it) - including the radial density distribution 
of stars around the binary. On the other hand, the flux 
of s tars into the loss cone is expected to peak around 
Th (jPerets fc Alexander! I2008D . For a given energy E 
and fixed M., J^ic{E,t) oc n{E,t)/T{E,t) where n{E,t) 
is the number of stars of energy E per unit (specific) 
energy, while t(E, t) is equal, in the spherical case, to 
Trix oc a^/p oc r^'"^!"^ (jBinnev fc Tremaini I2008D or, in 
the triaxial case, to the star's orbital precession time 

Turec(E,t) OC M*"^/^(< r) OC r'-^-^^/^, for r < few x rh 
()Yul I2002I ). Note that, for typical nucleus parameters, 
T'prec ^ Trix- As a rcsult, and since for more concen- 
trated nuclei (higher 7) n{Eh,t) is higher and T{E,t) is 
shorter, so more concentrated nuclei experience higher 
J^ic and higher s{t). 

The triaxiality also plays an important role in deter- 
mining the hardening rate. This is because of its role 
in supporting centrophilic orbits; in fact, at fixed n{E, t) 
only a fraction of the stars will be on centrophilic orbits, 
and this fraction is an increasing function of increasing 
triaxiality. As q decreases, the triaxiality of the remnant 
also decreases as shown in Figure ID This results pre- 
sumably in a weaker hardening rate since the fraction of 
stars on centrophilic orbits is smaller. This effect par- 
tially explains why the dependence of the hardening rate 
on M, is much weaker than expected based on the results 
obtained in Papers I and II for equal-mass mergers. 

We analyze the shape of the merger remnant by calcu- 
lating principal axis ratios and density profiles for galaxy 
merger models. Figure 0] shows the principle axis ra- 
tios of the merger remnant; these were defined as the 
axis ratios of a homogeneous ellipsoid with the same in- 
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Fig. 1. — Evolution of the separation between the SMBHs, in 
Af-body integrations of galaxy mergers for 7 = 0.5, 1.0, 1.5 and 1.75 
(top to bottom) according to Table [2] The arrows represent the 
time (T = 0) at which the two black holes become gravitationally 
bound for each model. Time and R are measured in A'^-body units. 



Fig. 2. — Evolution of SMBH binary semi-major axis, in N- 
body integrations of galaxy mergers for 7 = 0.5, 1.0, 1.5 and 1.75 
(top to bottom) according to Table [2] T measures the time since 
the instant when the SMBHs formed a bound pair. 
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Fig. 3. — Average hardening rates for 7 = 0.5, 1.0, 1.5 and 1.75 
according to Table [2] The average is measured between T = 60 
and T = 100 in A'^-body units. See text for details. 
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Fig. 4. — Ratio of intermediate to major (b/a, upper series of 
points) and minor to major (c/a, lower series of points) axes for 
7 = 0.5, 1.0, 1.5andl.75 according to Table[2] 

ertia tensor. The departures from spherical symmetry 
become more apparent for mergers with larger q. For 
q — 0.05 and 0.1, the merger-induced triaxiality becomes 
residual and the remnant results only a slightly flattened 
system. The dynamical effect of this transition seems to 
be abrupt, as the hardening rates shown in Figure [3] are 
essentially independent of mass ratio down to q 0.1. 
In such minor mergers, the secondary galaxy gets tidally 
disrupted within a few pericenter passages, and eventu- 
ally the naked black hole of satellite spirals in due to dy- 
namical friction. The evolution of SMBH binary in such 
a merger results very similar to those of binaries embed- 
ded in spherical galaxy models. Nevertheless in reality 
situations such as these are not very likely to happen 
as both galaxies would surely retain some triaxiality or 
flattening from a previous major merger. 

The density profiles of the newly formed galaxies as 
the result of the merger are shown in figure [S] As a 
result of the galaxy merger and the ejection of stars from 
the central galaxy region by the inspiraling black holes, 
the merger remnant nuclei have significantly shallower 
inner slopes than those of its progenitors. The "damage" 
caused by in-falling black holes increases with their mass. 
The values of 7t are calculated by fitting Dehnen's model 
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Fig. 5. — Density profiles of merger remnant for 7 = 0.5, 1.0, 1.5 
and 1.75 (top to bottom) according to Table |2] 
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to merger remnant right after the time when phase of 
rapid binary hardening ends. Table [3] presents the final 
inner slopes for merged galaxies. 

4. MASS DEFICITS 

Observations of nearby galaxies show that Sersic func- 
tions provide remarkably accurate fits for the major axis 
brightness surface profiles over the main bodies of bulges 
and early-type galaxies. This has been confirmed with in- 
creased accuracy and dynamic range over the last decade 
or so (Kormcndy ct al. 2009). Moreover, detailed state- 
of-the-art simulations of merging galaxies also gener- 
ate profiles which are consistent with Sersic functions 
(|Hopkins et al.l l2009al lH) . Even though there is not a 
complete theoretical understanding for the remarkable 
regularity of these profiles, their apparent generality and 
robustness led people to identify and interpret departures 
from Sersic fits and use them as a probe of the physics 
underlying the co-evolution of galaxies and their massive 
black holes. 

Cores tend to be found in giant ellipticals and are 
loosely defined as the central region in a bulge or early- 
type galaxy where the surface brightness deviates and it 
is below the values that would result from the extrapola- 
tion of the Sersic profile from the main body of the object 
down to its innermost region. Typically they are associ- 
ated to shallower cusps with 7 < 0.5 — 1. One possible 
explanation of a central core in a gas poor galaxy can 
be attributed to the ejection of stars by a hard SMBH 
binary. The destruction caused by the inspiraling mas- 
sive black hole in the center of galaxy can be measured 
in-term s of missing mass ca l led "m ass-deficit" . iGrahamI 
()2004[) and iFerrarese et al.l ()2006f ) measured the mass- 
deficits from the difference in fiux between the observed 
galaxies light profiles at the center and the inward extrap- 
olation from their outer Sersic profiles. Estimated mass- 
deficits from these studies yield Mdof '~ (1 — 2)M,, and 
values as high as ~ (4 — 5)M, were uncommon. More re- 
cent studies obtained larger values with an average value 

Mdef/M,) ^11 and an error in the mean of about 18% 

Kormendv eFaL| [2009^. 

Top panel of figure [6] shows the cumulative mass pro- 
file for model CI at different time of its evolution. The 
merger of a secondary galaxy with the primary causes 
some changes in the inner profile of the latter. This can 
be seen by comparing the profiles at times < = and 
t — 100, between which there is a noticeable drop in the 
stellar density around rh due to the merging of the two 
galaxies. Then, as the secondary is being tidally dis- 
rupted, after a few pericenter passages, the profile of the 
resulting merger remnant remains stable. Once, at time 
t ^ 202, the two black holes become bound and form a 
binary, the stars inside the loss cone start to be cleared 
by the binary and a rapid evolution of the mass profile 
ensues till t ~ 240. This drop in the central density is 
due to the slingshot ejection of stars by the binary. As 
a result, the the mass inside ~ few x rh drops rapidly 
hence changing the mass-profile there. Once this phase 
of rapid mass ejection ends, the density profile - notably 
inside rh - remains almost constant as the binary con- 
tinues its inspiral. This can be understood as follows. 
Inside rh, the gravitational potential is approximately 
spherical and thus the stars residing there can enter the 
loss cone only by diffusing in energy and angular mo- 
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Fig. 6. — Cumulative mass profile for one of our merger models 
at different time steps (top) and mass deficits normalized to M, 
for models C (bottom). Central density profiles at the beginning 
it = 0) and end {t = 280) are 7 = 1.5 and 7 ~ 1, respectively. The 
red arrow in the bottom panel of the figure shows the infiuence 
radius at the time of the formation of the SMBH binary. 

mentum space through two-body encounters. The cor- 
responding timescale is the relaxation time which scales 
as THx.j ^ ctVp ^ rT-3/2 (|Spitzeilll987b : since 7 < 3/2 
for all models at almost all times of interest, this implies 
that Tiix increases towards the center. The timescale for 
such stars to enter the loss cone is thus much longer than 
the typical precession time needed for centrophilic orbits 
of stars, resident at distances rh, to come sufficiently 
close to interact with the binary. As a result, the density 
profiles well inside rh remain essentially constant during 
the remainder of the inspiral so the damage impinged 
on the cusp by the SMBH inspiral following a merger is 
less severe than naive spherical model scenarios might let 
expect. The conclusion is that it is not very likely that 
a hole in the stellar distribution will ever be created by 
such mergers. 



We calculate the mass-deficits for each of our galaxy 
mergers as the difference between the stellar mass en- 
closed at a given radius at the time when the binary first 
becomes bound and at the end of the runs. The end of 
the run is a priori a somewhat arbitrary choice since not 
all runs end at exactly the same point in the binary's 
inspiral. However, when we compare the four cases (A2, 
A4, B3 and B4) for which we followed the binary up to 
final coalescence, we realize that most of the change in 
the density profiles happens at a relatively early stage in 
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TABLE 3 
Mass Deficit Analysis 







7f 


^final 


^ final 


^^-'dcf .1 


-''-'dcf,1.5 


^^^dcf ,2 


-''^dcf,2.5 


Af J r o 


Al 


0.19 


0.34 


6.4 X 10^* 


180 


0.51 


0.76 


0.84 


0.93 


0.72 


A2 


0.19 


0.26 


5.5 X 10~^ 


180 


0.50 


0.94 


1.20 


1.20 


1.05 


A2('PJ\f) 


0.19 


0.23 


TTiGV QG(i 


160 


0.48 


0.91 


1.11 


1.14 


1.01 


A '\ 


n 1 Q 


22 


^ Q V 1 


210 


u.uo 


1 21 


1 49 


1 40 


1 '\A 


A4 


0.22 


0.18 


9.7 X 10~^ 


220 


0.58 


0.76 


1.05 


1.50 


1.86 




0.22 


0.15 


1 1 LCI UCUj 


204 


0.59 


0.74 


1.06 


1.47 


1.81 


Bl 


155 


61 


2.9 X lO"^'^ 


140 


0.50 


1.18 


1.31 


1.38 


1.64 


B2 


14 




Q n V 1 C\^^ 


-LIU 


u.oo 


1 1 Q 


1 

i-.OO 


1 12 


i-.OO 


B3 


0.115 


0.51 


4.0 X 10"'^ 


115 


1.33 


1.73 


2.0 


2.13 


1.87 




0.115 


0.45 


1 1 LCI UCLL 


168 


1.69 


1.75 


2.10 


2.30 


1.91 


B4 


0.13 


0.48 


3.9 X 10~* 


230 


0.80 


2.0 


2.9 


3.50 


3.26 




0.13 


0.43 


merged 


189 


1.05 


2.05 


2.80 


3.34 


3.21 


CI 


0.066 


0.86 


1.25 X 10""' 


090 


0.71 


1.18 


1.42 


1.51 


1.29 


C2 


0.066 


0.80 


1.2 X 10-"' 


090 


0.96 


2.50 


2.88 


2.60 


2.40 


C3 


0.067 


0.77 


1.4 X 10""' 


080 


1.03 


1.52 


2.26 


2.40 


2.26 


C4 


0.069 


0.70 


1.02 X 10""' 


140 


1.10 


2.3 


3.9 


4.80 


5.30 


Dl 


0.045 


1.15 


8.3 X 10-5 


080 


0.91 


1.27 


1.82 


2.01 


1.89 


D2 


0.046 


1.06 


7.1 X 10-5 


100 


1.15 


1.85 


2.56 


3.04 


3.28 


D3 


0.046 


1.04 


7.0 X 10-5 


112 


1.39 


2.36 


3.14 


3.81 


3.97 


D4 


0.047 


1.01 


7.9 X 10-5 


130 


1.50 


2.84 


3.88 


4.96 


5.62 



Note. — In this table the columns from left to right represent: (1) Galaxy merger model, (2) the influence radius, (3)Thc inner density 
slope 7f for merger remnant, (4) the separation of the two black holes at the end of the run, (5) time at the end of the run in N-body units, 
(6—10) mass deficits M^^f „ in units of the mass of the binary Af, measured within n X rji ; rjj is calculated at the time of the formation 
of the SMBH binary 

the evolution of the binary, and this is mostly covered in 
almost all of the runs. The conclusion is that the derived 
mass deficit is quite insensitive to the exact final time of 
the runs, provided this occurs at a time when the binary's 
semi-major axis already reached a separation which is or- 
ders of magnitude shorter than the infiuence radius rh. 

Ours and the (several) observational definitions of mass 
deficit are obviously not equivalent, but they are surely 
somehow related provided the assumption that the mass 
is ejected by the binary holds true. We will not attempt 
to dwell into the intricacies of a detailed comparison, but 
will instead use the mass deficits obtained from our sim- 
ulations as qualitative indicators regarding the evolution 
of the galactic nuclei as they undergo minor /major dry 
mergers. 

The bottom panel of the Figure |6] shows the resulting 
mass-deficits measured out to various radii for model CI. 
Mass-deficits have a maximum value at around 2— 3rh for 
each model. We are showing the value of mass-deficit at 
different radii in table|3]for all the runs. First, we can see 
that the mass-deficits are clearly larger for more concen- 
trated (higher 7) models. They are also increasing with 
the mass ratio q. Note that what we are directly measur- 
ing from the runs is the effect of the inspiralling binary 
on the stellar distribution, not the total amount of mass 
ejected by the binary. The latter value should, by simple 
energetic arguments and given a fixed total binary mass, 
be on average the same regardless of the detailed proper- 
ties of the surrounding nucleus. The "damage" incurred 
by the stellar cusps in more concentrated nuclei, however, 
is expected to be greater as stars on loss cone orbits at 
small radii will represent a larger fraction of the (local) 
total stellar mass than those more outside. If galaxies un- 
dergo a series of mergers as suggested by the hierarchical 
galaxy formation scenario, then the mass deficits should 
add up because the time required to fill in the gap (re- 
moval of stars by inspiraling SMBH binary) is essentially 
of the order of relaxation time, which is several orders of 



magnitude g reater than a F lubble time for bright ellipti- 
cal galaxies (|Merrittll2006bl ). The follow up merger will 
redistribute the stars in phase space. The mass deficit 
following M dry mergers is M ■ {M^d / M,) so the esti- 
mated mass deficits obtained from our simulations imply 
that at least few (A/" =2 — 4) major me rgers are required 
to cre ate (Mdef/M,) ~ 11 reported bv iKormendv et all 
(I2OOI . 

It is also worth pointing out that the final density pro- 
file shown the top panel of Figure |6] has a central den- 
sity slope 7 1. This means that a single inspiral is 
not enough to turn a steep cusp (7 > 1.5) into a core 
(7 ^ 0.5). 

5. COALESCENCE TIMES FOR SMBH BINARIES 

In this section we want to derive an estimate for the co- 
alescence times of the SMBH bi naries in our simu l ations , 
similar to what has been done in lBerentzen et all (|2009t ). 
In the latter work it has been shown that the predictions 
of coalescence times for SMBH binaries in rotating, tri- 
axial galactic nuclei agrees well with the results of post- 
Newtonian A^-body simulations, up to order 2.5PN. To 
test the accuracy of our estimates made here, we repeat 
four of our A'^-body simulations of merging galaxies - this 
time including the post-Newtonian equations of motion 
for the SMBHs up to order S.5PJ\f (see Sec. [53]). 

In Fig. [21 one can see the rate of change of the binary's 
inverse semi-major axis during the hard binary phase is 
approximately independent of time, so we can write 



1 f da\ 

SNB = 2 I 7/7 ) ~ const., (8) 



a value which we can measure directly from each run. 

At some time in the simulation, the semi-major axis a 
will have fallen to such a small value that binary shrink- 
ing due to GW emission takes over. At such small separa- 
tions - depending on mass and eccentricity of the binary 
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- the gravitational wave emission becomes so efficient 
in extracting the angular momentum and energy from 
the binary that the latter becomes effectively decoupled 
from the rest of the stellar system and the SMBHs are 
led eventually to coalescence as if they were an isolated 
binary. 

The orbit-averaged expressions - including the lowest 
order 2.bVN dissipative terms - for the rates of change of 
a binary's semi-major axis, and ecc entricity due to GW 
emission are given bv lPetersI ([1961 : 



da 64 G^M,iM,2M. 

^dt^^^""Ta3c5(l-e2)V2 ^ 



1 



73 , 37 4 

— e -I e 

24 96 



,de, 304 G^M,iAI,2M, 



'dt' 



15 "a4c5(l 
121 
304' 



2)5/2 



1 



(9) 



(10) 



The corresponding hardening rate due to GW emission 
alone is thus given as: 



_ 64 G^M,iM.2M, 

. 73 2 37 4 

1 H H 

24 96 



(11) 



To calculate the coalescence time in our simulations we 
divide the e volution into two distinct regimes (see, e.g., 
iPreto et al.l ([2009)): (1) the classical regime, in which 
the hardening is driven by stellar-dynamical effects and 
(2) the relativistic regime in which the GW emission is 
dominant. We define this latter regime, starting from a 
time to when 

snb = Sow ■ (12) 

The semi-major axis ag at time can be calculated 
from Eq. [12] by assuming that both snb and eccentricity 
e{t) remain roughly constant^ during the stellar dynam- 
ical hardening phase, as supported by our simulations. 

We find that oq is typically smaller than afinai, the 
semi-major axis at the time ifinai at which our simulations 
end. Therefore, to usually exceeds ifinai and the time 
interval At between ifinai and to can be derived as 



At = 



'NB 



1 

ao 



1 



Ofinal 



(13) 



With this, the full time to coalescence, tcoai, in our 
(Newtonian) simulations takes into account (1) the time 
from when the two SMBHs become gravitationally bound 
(T) until they enter the GW regime to and (2) the bi- 
nary's lifetime icw until the final coalescence. The life- 
time tow of an isolated relativistic binary can be calcu- 
lated from Equations [S] and [TU] for a given semi-major 
axis oo and corresponding eccentricity eo, respectively 
(see Eq. 5.14 in Peters 1964). 

Note that if e(t) (slightly) increases as suggested by scattering 
experiment by Quinlan] 1119961 ) and Paper II then the full time to 
coalescence shall be smaller than our estimatedtcoal • Therefore, 
our results here should be interpreted as upper bounds to the true 
coalescence times. 



In order to compute the GW evolution rates one re- 
quires to adopt some physical units, e.g., of length and 
mass. In paper I, we equated the effective radii of our 
models to the effective radii of observed galaxies to fix 
the mass (Mgai) and the length (ro) units. In Paper II, 
we equate the influence radii of our models to that of the 
Galactic center and then extrapolate it to other galaxy 
masses by assuming the validity of the Af, — a relation 
and a stellar distribution p{r) cx r"^/'' at the center. 
Here to compute the coalescence times, we select three 
different elliptical galaxies (M32, M87, NGC4486A) to 
scale our models to physical units (Table |4]). We use the 
observed mass of the SMBH and its influence radius and 
compare it to the mass and influence radius of the pri- 
mary galaxy in our models to fix Afgai,p and ro.p. Model 
A's which have a very shallow slope are scaled with M87, 
model B's with NGC4486A and model D's are scaled 
with M32. The choice of galaxies for our models scal- 
ing is consistent with the fact that bright ellipticals have 
very shallow inner slopes and faint ellipticals have steep 
cusps. The estimated coalscence times for SMBH bina- 
ries in our simulations obtained through the proceduere 
described above are presented in Table [5l We find that 
the time interval spend in the Newtonian and in the rel- 
ativistic regime, respectively, are always comparable in 
almost all our simulations (see column 7 in Tab. [5]). This 
can be understood as follows: binaries with high eccen- 
tricities reach the relativistic regimes earlier than with 
low eccentricity and thus shortening the stellar dynam- 
ical phase. Also the higher eccentricities mean that the 
GW emission becomes more efficient which then short- 
ens the relativistic inspiral phase. The black hole of M32 
with a mass '-^ 3 x IO^AIq corresponds to possible sources 
detectable by LISA. The coalescence timescales for low 
mass black hole binaries are less than a Gyr (see table [5]) 
which suggests that prompt coalescence of binary SMBH 
detectable by LISA should be very common at high red- 
shifts. Even at the high mass end our models suggest 
the coalescence timescales ~ Gyr or even less depend- 
ing on the eccentricity of the binary. These timescales 
are short enough that SMBH binaries in these galax- 
ies should achieve full coalescence before a subsequent 
galaxy merger occurs. The coalescence times for SMBHs 
at low mass end are significantly (~ 10 x) shorter than 
the ones presented in paper I and more in line with those 
obtained in Paper II in which the scaling of lower-mass 
nuclei was set up by matching the total mass in stars 
to that observed in the inner ^ 2.5 pc of the Galac- 
tic center. The reason for this discrepancy is the weak 
dependence of effective radii with galaxy luminosity rela- 
tion adopted there, which was optimized for nuclei with 
> lO*Af0 black holes, but would lead to overestimate 
the influence radius of the Galactic center (or similarly 
compact nuclei) by an order of magnitude. Naturally, 
the coalescence times at high mass end of binary SMBHs 
are in nice agreement between both studies. 

5.1. Post-Newtonian Simulations 

In order to verify the accuracy of our estimate for the 
coalescence times given in Tabled we select the two cases 
A2 and A4 from Table [H] (plus two other cases, B3 and 
B4, not shown) and re-started these runs with exactly 
same initial conditions. The choice of these runs is moti- 
vated by their relative short coalescence times compared 
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TABLE 4 
Physical Scaling of our Models 



Model 


Galaxy 


M.{Mq) 


f'h(pc) 


T{Myr) 


L(kpc) 




speed of light (c) 


A 


M87 


3.6 X 10^ 


460 


1.9 


2.3 


7.2 X 10" 


257 


B 


N4486A 


1.3 X 10^ 


31 


1.4 


0.3 


2.6 X 10^ 


1550 


D 


M32 


3.1 X lO'^ 


3 


0.75 


0.12 


6.2 X 10** 


2011 



Note. — Columns from left to right; (1) Primary galaxy model, (2) Observed galaxy used to scale our model, (3) Observed mass of 
SMBH in the galaxy in column 2, (4) observed influence radius, (5)time unit, (6) length unit, (7) mass unit and (8) speed of light in model 
unit. 



TABLE 5 

Time to Gravitational Wave Coalescence 



Run 


"final 




^final 


eo 


ao (pc) 


to (Gyr) 


to /taw 


*coal (G 


Al 


6.4 X 10" 


-4 


9.10 


0.50 


3.5 X 10-^ 


1.30 


2.1 


1.89 


A2 


5.5 X 10" 


-4 


10.8 


0.98 


3.9 X 10" 


0.12 


1.1 


0.23 


A3 


5.9 X 10" 


-4 


10.3 


0.70 


6.9 X IQ-^ 


0.63 


1.2 


1.15 


A4 


9.7 X 10" 


-4 


9.60 


0.88 


1.6 X 10" 


0.30 


1.1 


0.57 


Bl 


2.9 X 10" 


-4 


23.3 


0.62 


7.4 X 10'^ 


2.10 


1.3 


3.70 


B2 


3.0 X 10- 


-4 


21.9 


0.98 


7.1 X 10-2 


0.24 


0.5 


0.77 


B3 


4.0 X 10" 


-4 


20.4 


0.95 


4.5 X 10-2 


0.39 


1.4 


0.66 


B4 


3.9 X 10" 


-4 


22.2 


0.96 


6.3 X 10-2 


0.27 


0.7 


0.64 


Dl 


9.2 X 10" 


-5 


75.7 


0.69 


2.1 X 10-^ 


0.48 


1.0 


0.98 


D2 


7.1 X 10" 


-5 


69.8 


0.66 


2.4 X 10-3 


0.43 


1.0 


0.82 


D3 


7.6 X 10" 


-5 


69.5 


0.61 


2.6 X 10-3 


0.41 


1.4 


0.70 


D4 


7.5 X 10" 


-5 


59.9 


0.60 


3.3 X 10-3 


0.38 


1.3 


0.67 



Note. — Columns from left to right; (1) Galaxy mergers model, (2) semi-major-axis (in model units) at the end of simulation tfinali (3) 
hardening rate (model units), (4) eccentricity at tflnali (5) semi- major axis of the binary at which the stellar dynamical hardening becomes 
equal to the hardening due to GWs, (6) Life time of binary in classical stellar dynamical hardening phase, (7) ratio between the time spend 
in classical regime and the time spend in GW regime and (8) Full time to coalescence of SMBHs. 



to other runs. 

We have implemented the relativistic effects to the 
MBH binary only by using the VJV equations of motion 
written in the inertial frame of binary center of mass 
including all the terms up to 'i.bVN order 



GM, 



[(l+^)ni2+Svi2] + 0(l/c8), (14) 



where n = r/r, the coefficients A and B are compli- 
cated ex pressions of the binary's relative separation and 
velocity ()Blanchetll2006D . The Post-Newtonian approxi- 
mation is a power series expansion in 1/c: the 0*'* order 
term corresponds to the dominant Newtonian accelera- 
tion. The IVN , 2VM and ZVM order terms are conser- 
vative and proportional to c^^, cr^ and c^^. The dissi- 
pative l.hVM and "i.hVM terms, which are proportional 
to and c~^, cause the loss of orbital energy and of an- 
gular momentum due to the gravitational wave radiation 
reaction. We treat the SMBHs as point particles and thus 
we neglect any spin-orbit or spin-spin coupling which in 
general is taken into account in the X.hVM (spin-orbit), 
IVM (spin-spin) and l.hVM (spin-orbit) terms in our 
common VM implementation. 

Figure [7] shows the evolution of the inverse semi-major 
axis 1/a and eccentricity e of the SMBH binary for run 
A2 with and without VM corrections. Hardening of the 
binary due to GWs emission starts to dominate at 1/a 
1000 which is slightly larger than the estimated 1/ao (all 
in model units). This is due to the slightly higher value of 
eccentricity reached during the run without VM which 
is then used to estimate ao. For the correct value of 
eccentricity obtained from run with VM , the estimated 
1 /ao matches accurately with the value where hardening 
due to GWs starts to dominate. The full coalescence 



time for the binary in our VM simulations is 0.25 Gyr 
which is very close to to estimated tcoai ^ 0.23 Gyr. 
The inverse semi-major axis and eccentricity evolution 
for model A4 are shown in figure [S] If we look at the 
evolution of the inverse semi-major axis for this run then 
we see that hardening due to GW becomes important 
at estimated 1/ao, as the average eccentricity evolution 
matches almost perfectly between the VM and mm-VM 
runs. 

As in the (isolated) models of iBerentzen et all ()2009D . 
we find very good agreement between our estimated co- 
alescence times and those found using full VM terms in 
our merger runs. We should note that the coalescence 
times for Models D are all in the upper range of the co- 
alescence times estimated in Paper II for lO^M© LISA 
binaries. This can be easily understood if we note that 
in the previous paper, the nuclei were all modeled with 
7=1, whereas here we choose 7 = 7/4 more in line with 
the expected properties of compact nuclei. We can see 
both from Figure [T] and Table [5] that the binary eccen- 
tricities tend to be lower (e ~ 0.6 — 0.65) for 7 = 7/4, 
and indeed in Papers I & II, the binaries reached very 
often values e > 0.9 (in accordance to values obtained 
in this paper for 7 = 1). The dependence on eccen- 
tricity of the coalescence time under GW emission is 
,GW ^ (1 — e^)'^/^ could easily account for a de- 
crease of an order of magnitude or so when the binaries 
are very eccentric. It will be accordingly very important 
to further investigate the dependence of the eccentricity 
evolution under different values of 7 and g, and this is the 
subject of a forthcoming paper (jPreto et aLllin prep.D . 

6. SUMMARY & CONCLUSIONS 

Direct A^-body simulation of mergers of spherically 
symmetric galaxies with different mass ratios were per- 
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Fig. 7. — Evolution of the inverse semi-major axis 1/a (top) and 
eccentricity e (bottom) for model A2, with and without VM terms. 
The horizontal lines represent the estimated semi-major axis of the 
SMBH binary for which the stellar dynamical hardening becomes 
equal to the VM hardening derived from the run without VM (ao) 
and with VM (ao — VM). See main text for further details. 

formed to investigate the evolution of binary SMBHs 
from the onset of merger through the stellar hardening 
phase until the eventual relativistic coalescence. The 
merging galaxies have different initial density profiles 
varying from shallow (7 = 1/2,1) to steep cusps (7 = 
3/2,7/4) - thus covering the range of stellar distribu- 
tion typically observed in the centers of bulges or early- 
type galaxies. In Papers I and II, we have shown that 
merger-induced triaxiality could support a purely stellar 
dynamical solution to the FPP in the case of equal-mass 
mergers. In order to assess the prospects for such a solu- 
tion in the more general unequal-mass case, we measured 
both the hardening rate and merger-induced triaxiality 
for mass ratios in the range q g [0.05, 1]. The merger- 
induced triaxiality found in Papers I and II for equal- 
mass mergers is still present in unequal (g < 1) mergers, 
albeit it becomes weaker as q decreases. Minor merg- 
ers with q < 0.05 of spherical progenitor galaxies leave 
the primary almost unperturbed, so the subsequent bi- 
nary evolution follows suit in an almost spherical back- 
ground nucleus. The classical FPP could well show up 
in such cases. This transition seems to be abrupt, as the 
measured hardening rates are essentially independent of 
the mass ratio q until it suddenly declines at a value 
somewhere between q = 0.05 and 0.1 - thus indicating 
that only a very modest triaxiality is needed for driving 



. 8 



0.4 



0.2 




250 



Fig. 8. — Evolution of inverse semi-major axis 1/a (top) and 
eccentricity e (bottom) for model A4 with and without VAf terms, 
ao is the semi-major axis of the SMBH binary for which stellar 
dynamical hardening is equal to VAf hardening. 



the binary inspiral at rates consistently higher than in a 
spherical nucleus. It is also not surprising that we find 
that the hardening rate increases substantially for high 
7, as more concentrated nuclei will have bigger amounts 
of centrophilic orbits. 

The solution to the FPP therefore hinges strongly on 
the expected triaxiality - or, more precisely, the amount 
of centrophilic orbits supported by it - at the center 
of galaxy merger remnants. From this perspective, our 
models may represent worst case scenarios^ since they 
were derived from the merger of spherically symmetric 
progenitor galaxies, which is indeed somewhat unlikely 
as most observed ellipticals are significantly fiattened 
and a good fra ction of them are at least mildly triaxial 
(IKormendv fc B ender 1996; Emscllem ct al. 2007). How 
the shape of the progenitor galaxies affects the shape 
of the merger remnant and the resulting SMBH hard- 
ening rates is yet unclear. Nevertheless it is important 
to pursue increasingly realistic stellar dynamical mod- 
els of galactic nuclei in order to refine our understand- 
ing of binary SMBH evolution tracks and associated 
timescales. A detailed study of the A^-independence and 
merger-induced triaxiality in unequal-mass mergers - in- 
cluding non-spherical progenitors resulting fr om previous 
merger e vents - is being currently pursued (jPreto et al.l 
lin prep. I). 

We have produced estimates of coalescence times us- 
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ing a simplified prescription for the late relativistic phase 
of the inspiral - i.e. adopting the Peter's (1964) equa- 
tions for the evolution of an isolated binary and ignoring 
the late stellar-driven eccentricity evolution. We assessed 
the accuracy of these approximations by including Post- 
Newtonian terms (up to 3. SPA/" order) to the equations 
of motion of the binary in a few representative cases of 
our sample. At least in these few selected cases, the 
agreement is remarkably good, and this is mainly due to 
the fact that the eccentricity evolution is almost iden- 
tical in the VJV and non-'PA/' runs. One should add a 
little note of caution though. The eccentricity evolu- 
tion may be a quite sensitive function of the properties 
of the surrounding stellar cluster (slope 7, mass ratio q, 
amount of net rotation, etc) and on the initial conditions 
of the binary. Previous studies indeed indicate that one 
could expect significant evolution of the eccentricity dur- 
ing the har dening phase , even though not always of the 
same sign (|Sesanall201Cll : iPreto et"aLl 120111 : iSesana et al.l 
|2011| ). The coalescence times resulting from our calcu- 
lations are consistently shorter than a Hubble time at a 
given redshift (up to 2: 2 — 3 for all masses and up 
to z ~ 6 — 10 for those in the LISA mass range). Co- 
alescence times for binaries of ^ IO^Mq are all shorter 
than 1 Gyr; for the most massive ~ 10 Mq they range 
between hundreds of millions of years (highly eccentric 
ones) to ~ 1 — 2 Gyr (less eccentric ones). SMBH bi- 
naries are therefore very promising sources of GWs for 
LISA both at low and high redshift. The coalescence 
times obtained here - especially those from models D 
- may be effectively an upper bound if the mild eccen- 
tricity growth observed in these runs is not generic (we 
have only four runs). In Paper II, we used models with 
7 = 1 to compute the evolution of binaries in the LISA 
mass range and we obtained faster coalescences; this was 
partly because the binaries were significantly more eccen- 
tric. Therefore the detailed properties of the distribution 
of coalescence times hinges in part on the dependence of 
the eccentricity evolution on the properties of the sur- 
rounding cluster. The question of whether eccentricity 
grows, decays or remains approximately constant dur- 
ing the hardening phase remains still uncerta in and we 
are currently pursuing it (jPreto et aLllin prepi ). Interest- 
ingly, the duration of how long the binaries remain in the 
stellar dynamical hardening phase and in the relativistic 
one are roughly equal by an eccentricity regulated pro- 
cess. This means that SMBH binaries on highly eccentric 
orbits have less probability to be detected in the classi- 
cal regime but at the same time are sources of strong 
GW signals to be detected. On the other hand, binaries 
with low eccentricities remain longer in the Newtonian 
regime and therefore are more likely to be observed in 
this phase, but may only be luminous GW sources in the 
final phase of coalescence. 

With our results we moved another step towards a con- 
sistent stellar dynamical solution to the FPP, and thus of 
providing a solid dynamical substantiation to cosmolog- 
ical scenarios where prompt coalesce nces are the norm 
during SMBH- galaxy co-evolution (jSesana et al.l 120071 : 
IVolonteril[2010l ). If our results hold true in real galaxies, 
the bottleneck to SMBH coalescences - if any - is likely to 
be associated to the long timescales needed for SMBHs 
to become a bound pair in galaxy mergers of unequal- 



mass -_es£eciall5;_Jii case g ^5 0-1 s-nd gas fractions are 
low (iCallegari eralll2011[ ). These results are promising 
for SMBH binaries being abundant LISA sources at high 
redshift. 

What about the role of gas in hardening the SMBH 
binaries at sub-parsec scales? Numerical simulations of 
merging disk galaxies with a substantial gas- fraction have 
shown that gas is driven towards the central region of 
the merger remnant where it accumulates. As a con- 
sequence, the remnant becomes more axisymmetric (i.e. 
losing triaxiality) and the orbital structure changes no- 
ticeably as compared to the purely stellar case (Naab, 
Jesseit & Burkert 2006). Thus we can expect that the 
fraction of centrophilic orbit is affected in the same man- 
ner by the presence of the gas. In which direction it 
affects the hardening rate is yet an unresolved issue. 

The gas which is funneled to the center of the remnant 
may undergo strong starbursts which may lead to the 
formation of a steep stellar density cusp - even before the 
binary has formed - as shown by Callegari et al. (2009), 
which then in turn increases the binary's hardening rate 
as shown in this work. 

The gas could certainly also assist the inspiral if it 
is cold enough to settle into a circumbinary thin disk 
(|Armitage fc NataraianI [2005D . even though the uncer- 
tainties associated with its long-term dynamical behav- 
ior are potentially more severe than those related to stel- 
lar dynamics. For instance, gaseous disks are suscepti- 
ble to fragmenting and forming stars. In order for the 
disk to be at worst marginally stable against fragmenta- 
tion at every radius, it s total mass is constrained to be 
Mdigk ~ (0.1-0.2)M. (jCuadra et al.ll2009l:lLodato et all 
|2009| ). This is not likely to be enough to drive the binary 
to coalescence. 

The "mass deficits" induced to the centers of galaxies 
with no gas by inspiraling SMBH binaries are found to 
be of order ~ (1 — 5)M,, depending on the slope 7 of 
their central stellar distribution and on their mass ra- 
tio q. This is consistent with the picture that cores in 
giant ellipticals are scoured by SMBH binary inspirals 
during successive merger events - assuming only that 
relaxation times are long enough in those systems for 
the results from different mergers to be cumulative. For 
galaxies in the mass range of the Milky Way the picture 
is basically distinct. Since the amount of mass deple- 
tion in a single merger event only partially destroys a 
steep inner cusp, it results that even if a major merger 
were to have happened in the Galactic center at redshift 
2: = 1 or larger, it would have been enough time to re- 
grow a mass segregated Bahcall & Wolf cusp since then 
(jPreto fc Amaro-Seoanil2010f ). 
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